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Bosonic van der Waals clusters of sizes three, four and five are studied by diffusion quantum 
Monte-Carlo techniques. In particular we study the unbinding transition, the ultra-quantum limit 
where the ground state ceases to exist as a bound state. We discuss the quality of trial wave functions 
used in the calculations, the critical behavior in the vicinity of the unbinding transition, and simple 
improvements of the diffusion Monte Carlo algorithm. 
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I. INTRODUCTION 



In a previous paper0, a form of variational trial wave function was developed and applied to atomic, few-body 
5—1 ' systems consisting of five or less atoms of Ar and Ne interacting via a Lennard-Jones potential. In addition, we tested 
the trial functions for a hypothetical, light atom resembling Ne but with only half its mass. We did not study atoms 
such as ^He with larger de Boer parameters, i.e., systems in which the zero point energy plays a more important 
, role relative to the potential energy. Studying such systems is the main purpose of the present paper. In fact, we 
study clusters in the ultra-quantum unbinding limit, in which the zero-point energy destroys the bound ground state. 
\ Simple arguments applied to this unbinding transition predict the way in which the energy vanishes as the de Boer 
^ • parameter approaches its critical value, and the nature of the divergence of the size of the clusters. Our numerical 
' } ' results are in agreement with these predictions. 

, As the de Boer parameter increases, the quality of the wave functions decreases and, whereas in Ref. |^ there 
• was no need to go beyond variational Monte Carlo, we rely in this paper on diffusion Monte Carlo to improve the 
^ variational estimates. Since the diffusion Monte Carlo algorithm is based on a short-time expansion of the imaginary- 
\^ time evolution operator exp(— rTi), where Ti is the Hamiltonian, the algorithm is exact only in the limit r — > 0. In 
0^ ' practice, since the computations are done at finite r, and the computer time required to obtain a given statistical 
i accuracy increases as t~^, it pays to design the diffusion Monte Carlo algorithm to have a small time-step error. For 
^ this purpose we use a simplified version of the improved diffusion Monte Carlo algorithm introduced by Umrigar et 
' ala We use the trial wave functions of Ref. |l| to obtain both variational and diffusion Monte Carlo estimates for the 
^ ' ground state energy and other expectation values. 

D , The layout of this paper is as follows. In Section || we review diffusion Monte Carlo and the modifications we 
' made to the algorithm given in Ref. ^ to make it applicable to Lennard-Jones bosons. Trial functions are optimized 
^ ' by minimizing the variance of the local energy. In the immediate vicinity of the unbinding transition this method 
becomes unstable, a problem that can be solved straightforwardly, as discussed in Section We found that graphical 
methods were helpful to identify configuration space regions that dominate the variance of the local energy. Some of 
the details are discussed also in Section O. However, this part of our presentation is incomplete in the sense that our 



' techniques rely in part on color graphics, which are not included in this paper. For additional details we refer to Ref. ||. 
Section III contains numerical results: quantitative evaluation of the improvements obtained by the modified diffusion 
Monte Carlo algorithm; estimates of ground state energies; and numerical corroboration of the critical behavior of 
the unbinding transition, as predicted below [c/. Eqs. (jig)]. 



II. METHODS 



We consider a clusters of N bosonic Lennard-Jones atoms. The Lennard-Jones pair potential v{r) ~ Ae[{r / a)^^"^ — 
(r/cr)"^"] in reduced units takes the form v{r) — r^^^ — 2r~^, and the only independent parameter in the Schrodinger 
equation is the reduced inverse mass m~^, a quantity proportional to the square of the de Boer parameter, h/ay/me, 
a dimensionless quantity measuring the importance of quantum mechanical effects. 

A cluster configuration is given by the cartesian coordinates of the atoms, which form a 3A^-dimensional vector 
R = (ri, . . . ,rjv), where is a 3-dimensional vector specifying the coordinates of atom i. The total potential energy 
of the cluster is denoted by V(R). 

We briefiy review the diffusion Monte Carlo algorithm, implemented, as usual, with importance sampling for which 
an optimized trial function ?/'t(R-) is introduced. The variational energy of this state is written as Et- For a given trial 
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wave function '0t(R), one introduces a distribution, /(R, i) = ?/jt(R)V'(Ri i)- Since ipCRjt) satisfies the Schrodinger 
equation in imaginary time, /(R, t) can be shownuu to be a solution of the equation 

- VV(R, t) + ^y- [V(R)/(R, t)] - 5(R)/(R, t) = (1) 

where the velocity V, usually called the quantum force, is given by 

V(R) - (vi, . . . , v^) - -^-^ , (2) 

and the coefficient of the source term is defined as 

5(R) = Et- £(R), (3) 

which in turn is defined in terms of the local energy 

7^^t(R) _ 1 V^^t(R) , . , . 

" - '2^ ^t(R) + ^ ^ 

Note that the energy in the Schrodinger equation was shifted so that S vanishes if the trial state is an exact energy 
eigenstate. 

The simplest implementation of the diffusion Monte Carlo algorithm is based on a short-time approximation of 
G(R',R, r) = (R| exp(— T7i)|R'). This takes the form of a product of the Green hmctions of each of the three 
operators on the left-hand side of Eq. (^: 

^(R', R, rm) « (27rT)"^ j dR'^e' ""''^^"'^ ^[R" - R - V(R)T]e-^[^(^')+^(^)l™. (5) 

The Monte Carlo incarnation of the above expression consists of a deterministic drift of an initial configuration R 
to a new configuration R" followed by a random diffusion to a final configuration R'. Finally, there is a reweighting 
based on the initial and final configurations. According to Eq. (|^), both drift and diffusion modify the coordinates of 
each atom in the configuration simultaneously. Alternatively, one can break up the operators on the left-hand side 
of Eq. (P into single particle operators. Correspondingly, one can write a short-time Green function as a product 
of factors associated with drift and diffusion of individual atoms, e.g., in the order defined by the numbering of the 
atoms i = 1, . . . ,N . In our computations we used this alternative approximation, while we kept the same exponential 
growth-decay factor as in Eq. rather than including a reweighting factor for each single-atom update. 

The imaginary-time evolution operator does not uniquely define a short-time expansion. The corresponding freedom 
can be exploited to extend the range in r over which the approximate short-time Green function agrees with the exact 
expression to some given accuracya. In other words, the time-step error of the diffusion Monte Carlo algorithm can be 
reduced by adapting the algorithm so that it can deal more accurately with singular regions of configuration space. In 
this way, one can make simple algorithmic changes of essentially zero computational cost to improve the efficiency of 
the algorithm dramatically. Indeed, this was the guiding principle in the design of the diffusion Monte Carlo algorithm 
described in Ref. ^ 

When one is dealing with atoms or molecules in which the "elementary" particles are the electrons and nuclei, 
there are numerous sources of singular behavior: electron-electron and electron-nucleus collisions and nodes in the 
trial function. In the bosonic case, the ground state has no nodes and the only singularities are due to inter-atomic 
collisions. For example, the exact Green function of the second term, i.e., the drift term in Eq. (|^) is (5[R" — R(t)], 
where R(i) is the position at time t obtained by exact integration of the velocity V subject to the initial condition 
that R(t) = R at t = 0. This exact expression reduces to the approximation (5[R" — R — V(R)r] if V is assumed 
constant during the time r in Eq. (^, but this assumption fails when two atoms collide, as can be seen as follows. 

One of the boundary conditions one would like to impose on the trial wave function is that the local energy remain 
finite when the distance between two atoms vanishes, but unfortunately, we have not been able to meet this goal. 
In particular, suppose that the distance r between two atoms vanishes while all other distances remain finite and 
non-zero. By choosing a trial wave function that behaves as exp(— -y/ro/Sr^l, one can satisfy the condition that for 
r ^ the local energy does not diverge as strongly as the potential energjtJ, i.e., as r^^^, but only as r^^. In this 
case, Eq. (||) implies that the velocity diverges as r~^\/rn. This divergence implies that for sufficiently small r the 
approximation R(t) w R -|- V(R)t becomes a poor one, since it is obtained from the assumption that the velocity is 
constant during the time r. To improve the approximation following Ref. |2| we integrate the speed as given by the 
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differential equation v = r ^ ^pm for the two-particle problem and express the resulting average speed in terms of the 
initial speed. Applied to the drift of atom i this yields 



_ -l+\/l + 7(;^)^jT 

vii = —_ Vi. (6) 

V m / ^ 

If ViT, the one-particle drift for atom i, is replaced by ViiT, the original expression is reproduced for small velocities, 
while for large velocities the magnitude of the drift is reduced to (jTyfm)^!'^ . 

The problem of the short-range singularity is most pronounced for light particles. For heavy particles, on the 
other hand, the wave function is strongly peaked close to the classical equilibrium position Rq, and in this case the 
approximation of constant velocity can be improved too. Suppose we assume a Gaussian approximation for the wave 
function 

V't(R) oce-^(^-^°)', (7) 

where Rq represents the classical configuration of minimum energy. In this approximation, the velocity is always 
directed towards Rq and vanishes at that point. To compute approximately the drift of atom i for a trial function of 
this form, one can express A in terms of the local kinetic energy of particle i and its velocity, given by Eq. (H) , as 

^-6^"^' ^t(R) ^' ^ ^ 

an expression containing only quantities that have to be computed anyway. As in the case of the diverging velocity, 
one can integrate the velocity exactly and express the result in terms of an average velocity 

1 _ £-2^^ 

V2,: = (9) 

The diffusion Monte Carlo algorithm cannot be expected to be efficient unless values of r are chosen so that the 
order of magnitude of a typical drift or diffusion step is comparable to the width of the region in which the wave 
function is appreciable. In the classical, large m limit. At k, mr ~ 1. This implies that in practice expression (j^) 
for the mean drift velocity will not differ dramatically form the one given by the original expression (^). In fact, we 
found that the impact of this modification of the algorithm on the reduction of the time-step error and increase of 
efhciency was insignificant. 

We made no attempt to construct a sophisticated scheme to interpolate between the two approximations as given 
in Eqs. and (^). Instead, we simply used the average velocity in all our computations: 

= min(vH, V2i), (10) 

where value of the function min is the vectors with the smallest magnitudes. 

Finally, in order to guarantee that the diffusion Monte Carlo algorithm produces the exacl Green function in the 
ideal case that the trial function is the exact ground state we include an accept-reject stepd. Once all atoms have 
drifted and diffused, a new configuration R' has been generated. This configuration is accepted with probability 

. JVT(R')l'G'(R.R',r) , , , 

p ^ min {^-^ — ^' , ^ ' 1 . (11) 

'|^T(R)PG(R',R,r)' 

If the new configuration R' is rejected, the previous configuration R is kept. We note that the accept-reject step 
requires for its implementation the introduction of an effective time step TcS in some parts of the algorithm. For more 
details see the electronic structure algorithm in Ref. H, which has to be simplified in obvious ways to be applicable to 
the current, atomic system. 

For a given amount of computer time, the statistical accuracy of the diffusion Monte Carlo computations can be 
increased by using optimized trial functions. We used trial wave functions described in Ref. |l|. They were optimized 
by minimization of , the variance of the local energy, but we found that this optimization procedure was not stable 
close to the unbinding transition. This instability can be understood as follows. 

The variance of the local energy cannot be evaluated exactly, or even numerically exactly, for an arbitrary trial 
state, since this would require a 3n-dimensional integration. Instead, one uses a Monte Carlo approach in which a few 
thousand states are sampled from the square of the trial wave function defined by an initial guess for the parameters 
to be optimized. Then, one changes the parameters and estimates the variance by reweighting configurations with the 
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appropriate ratio of the current probability density and the one from which the sample was drawn originally. This 
can be done efficiently as long as the the two wave functions have sufficient overlap. Once this condition is no longer 
satisfied, one generates a new sample from the current distribution and iterates until the process converges. 

The energy of clusters goes to zero when the de Boer parameter approaches a critical value, and at the same 
time these clusters grow in size. Under these circumstances, Monte Carlo samples of fixed size tend to consist 
of configurations predominantly sampled from the tail of the wave function. During the optimization of the wave 
function, the local energy tends to become the same, physically incorrect constant for these configurations, and as 
a consequence the variance of the local energy as estimated from a sample of a fixed size can be reduced artificially 
by choosing an energy even closer to zero. Of course, the true variance of the trial wave function increases in this 
process, but for a sample of fixed size this goes undetected. We found that the solution to this instability is quite 
simple: rather than minimizing one minimizes x^/£'^. 

Computations of the ground state energy by variational or diffusion Monte Carlo satisfy a zero-variance principle, 
i.e., in the limit where the trial function is the exact ground state eigenfunction, the statistical error vanishes. As a 
practical consequence, computations can be made considerably more efficient by a good choice of trial functions. In 
the design of trial wave functions described in detail in Ref. |^, we followed the same procedure as in Refs. Hj^: the 
trial functions satisfy boundary conditions associated with (a) the collision of two atoms and (b) having one atom go 
off to infinity. The most likely configurations, which involve intermediate distaaces and require most of the variational 
freedom of the trial wave functions, are described by many-body polynomials. 

In the process of improving the quality of wave functions it is important to know what region of configuration space 
contributes most to the variance of the local energy. For instance, it is useful to know if the quality of the wave 
function is limited by poorly satisfied boundary conditions or whether the quality can be improved by adding more 
variational parameters. Another possibility is that the wave function has too much variational freedom relative to 
the sample over which it is optimized. This might lead to unphysical peaks in the wave function, which might only 
show up in the variance of the local energy obtained from production runs, which sample a much larger number of 
configurations than the number present in the sample used for the optimization of the trial function. 

To help answer such questions, we made density plots of the local error, the deviation of the local energy from its 
average. As an illustration, we discuss the case of five-atom clusters. In fact, we used superimposed color density 
plots of both the wave function and the local error, which contain more information than can be reproduced by the 
grey-scale plots reproduced in this paper. We refer the reader to Ref. ^ for the color graphics. 

Obviously, the fact that the ground state wave function depends on 3N — 6 independent coordinate variables, 
seriously limits any graphical approach. For the five-atom clusters we found the following planar cut through config- 
uration space informative: four atoms were fixed at the vertices of a regular tetrahedron, while the fifth particle was 
located in a plane that contains two of these vertices and bisects the edge connecting the two remaining atoms. 

Figs. and || strongly suggest which regions of configuration space contribute most to for case N = 5. For 
the interpretation of the density plots the following convention should be used: zero intensity (white) corresponds 
to a minimum, while full intensity (black) corresponds to a maximum of the plotted function. Fig. [l] represents the 
density plot of the weighted "local error", defined as \{£ — -Et)V't|, [c/. Eq. (^] as a function of the position of 
the fifth, wandering atom. Note that the quantity to be minimized in the optimization of the trial functions is 
the configurational integral of the square of this quantity, apart from a normalization constant. Fig. ^ shows the 
dependence of on the position of the fifth particle while the other atoms are fixed in the geometry discussed 
above. 

The conclusion we draw from the density plots is that the trial function fails particularly in regions where more than 
two atoms collide and we see that of these the the local error is largest whenever four atoms are close. Unfortunately, so 
far we have not been able to find trial functions without this problem, i.e., trial functions without the r^^^ divergence 
of the local energy which occurs when more than two atoms collide, and without the r^^ divergence of two particles 
in the even distant presence of a third. 

III. RESULTS 

A. Groundstate energy and time-step error 

We present results for the time-step error, discussed in Section || and compare the version of the diffusion Monte 
Carlo algorithm summarized above with a simple version of the algorithm in which (a) the velocity is treated as a 
constant for the integration of the short-time drift Green function; (b) each move is unconditionally accepted rather 
than the result of an accept-reject step, so that Tcff — r. 
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FIG. 1. 

Density plot of the "local error" in the geometry described in the text. The two dots in the lower right-hand corner arc the two in-plane 
vertices of the tetrahedron; the one in the upper left corner is the projection of the two out-of-plane vertices. The length of the tetrahedron 
edges is 1.3 and = 0.16. The darker the region, the more it contributes to x^- Note that the dark region in the lower right is a cut 
through the banana-shaped dark region in the upper left. The regions of the largest local error are the two symmetrically located regions 
where the wandering atom is close to three others. White lines are are cuts through nodal surfaces of the local error and have no physical 
significance. 
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FIG. 2. Density plot of |V't|^ as a function of the position of the fifth particle in the geometry described in the text and the figure 
caption of Fig. ^. Dark regions correspond to high probability density. Note that the regions with the largest local error are contained in 
the region where the wave function becomes small because of repulsive potential at short range. 
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FIG. 3. Estimate of the ground state energy Eo as a function of time step for Nes 

We compared the simple and improved diffusion Monte Carlo algorithms for clusters of Ar, Ne and hypothetical 
"^-Ne" atoms with sizes in the range = 3, 4 and 5. Figs. ||, and ^ display plots obtained for estimates of the 
ground state energy Eq for Ne and "i-Ne" clusters with N = 5; the behavior for the smaller systems is analogous. 
As expected, the reduction of the time-step error is greatest for the lighter atoms. For Ar no significant improvement 
of the algorithm was achieved, and only an approximate reversal of the sign of the time-step error occurred. 

Table | displays estimates of the ground state energy obtained by variational Monte Carlco and the improved 
diffusion Monte Carlo algorithm. In addition, the table contains information pertaining to the magnitude of the bias 
of the variational estimate due to the fact that Et is an upper bound of the true ground state energy. The tightness 
of this bound is determined by the quality of the trial wave function. More in detail, the variance of the local energy 
is defined by 



2_ {jjTlin ^ Eofl^T) 

(V'tIV't) 

and the following inequality holds (see Ref. |l] for details and references): 



< £:t - £^0 < 



El — Eo' 



(12) 



(13) 



where Ei is the energy of the first, totally symmetric excited state. To estimate the number of correct digits in the 
variational estimate of the ground state we use the following quantity: 



Q' = - log 



X 
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(El — Et)\E't\ 



(14) 



It is also of interest to know how tight a bound the right-hand side of inequality (|l^ ) provides. This is measured by 
the following ratio 



R = 



{El — Et){Et — Eq) 

The results are shown in Table |[ Quite remarkably, the bound given in Eq. ([l3| ) is very tight. 



(15) 
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FIG. 4. Estimate of the ground state energy Eo as a function of time step for ^-Nes 



TABLE I. Diffusion Monte Carlo estimates of the ground state energies Eq for noble gases Ar and Ne, and hypothetical 
"i-Ne", compared with variational Monte Carlo estimates taken from Ref. Standard errors in the last digits are given in 
parentheses. 
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FIG. 5. Fitting curves for diffusion Monte Carlo estimates of the ground state energy for clusters of sizes A'^ = 3, 4 and 5. 



A severe test for the accuracy of a cluster trial function is its performance in the strong quantum limit, i.e., for 
large values of the de Boer parameter. In particular, we discuss results in the vicinity of the unbinding transition, 
where the cluster ceases to possess a bound state. 

The ground state of the ''He dimer is believed to be a (weakly) bound state. Since the ground state energy 
presumably decreases with cluster size and other systems have smaller de Boer parameters, we can safely assume that 
the unbinding transition for boson clusters is inaccessible experimentally. However, the transition does occur at finite 
cluster size for ^He, and it makes sense to use the boson case as a simpler test case for the trial functions. 

A second issue of theoretical interest is the behavior of energy and cluster size as a function of the dc Boer parameter 
in the vicinity of the unbinding transition. 'Phis transition plays the role of a critical point and, in fact, has many 



The following critical behavior is expected for the ground state energy iJg and the average size (r) (as defined 
below) of the cluster 



for m I 7Tic where Am = m — ra^ with the critical value of the dimensionless mass. 

The critical behavior given in Eq. (|6|) can be made plausible as follows. For the simple case of a dimer one can 
show this directly, and the mathematical mechanism that yields Eqs. (|l^) is the following. Two scattering states 
wave functions forming a complex conjugate pair merge at zero momentum to produce two states with "complex 
momentum" : a physically acceptable bound state and a state with unacceptable behavior at infinity. This mechanism 
is probably not limited to the dimer, and therefore it is quite plausible that Eqs. (^6|) apply in general to clusters of 
any finite size. 

On the other hand, in the m — > cxd limit (vanishing de Boer parameter), the harmonic approximation predicts 



is expected to be linear in the de Boer parameter in both the classical and extreme quantum limits. Indeed, figure o 
displays this remarkably dull behavior for the square root of the normalized energy as a function of the de Boer 
parameter over the whole range. To display the critical behavior of the energy in more detail. Figure ^ shows a double 
logarithmic plot of E versus m — m^. 

Numerical values for average distance and the gyration radius were obtained using the approximation 



B. Unbinding transition 




-Bo - (Am)2, 
(r) - (Am)-i 



(16) 
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FIG. 6. Double-logarithmic plot of Eo for A'' = 3,4 and 5 versus m — trie. 



(r) w 2(r)^T,o - (r)^T,i>T> (18) 

where the first term on the right denotes the mixed expectation value obtained by diffusion Monte Carlo, while the 
second term is the variational expectation value of the cluster radius in the trial state. 

Figs. 1^ and |^ are plots of the approximate values of the average cluster size, as measured by the average inter-particle 
distance and the gyration radius, versus inverse m — rric. The behavior displayed in the graphs is consistent with 
the scaling law given in Eq. (|l6|). It should be noted that there are apparent irregularities in the data points. These 
can be traced to irregularities in the quality of the wave functions, which are a result of incomplete optimization. 
As far as the energy is concerned, the quality of the trial wave functions only affects the statistical accuracy of the 
estimates, but, as can be seen from Eq. |l^, imperfections of the optimized trial wave functions result in true errors of 
expectation values of quantities that do not commute with the Hamiltonian, such as the size of the clusters. 
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FIG. 7. Double-logaxithmic plot of the approximate the average size of a cluster (r) defined as the average inter-atomic 
distance versus m — rric for A?^ = 3, 4 and 5. 
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FIG. 8. Double-logarithmic plot of the approximate values of the average size of a cluster (r) defined as the gyration radius 
versus m — rric for iV = 3, 4 and 5. 
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